% section 6.4.1 'robustness' / Different starting states 
     % this file shows the differences in CE for GHHW and GLTHI 
     % for many draws of starting probabilities

    clear;
    addpath functions
    close all
    rng(1234567,'twister')

    load('./results/model_parameters_main.mat');
    load('./results/capHtmod_c.mat');
    % file sampler_probs_norm_c.csv used in the paper not provided 
    % in replication package due to its size 
    replicate = 0;

    if replicate==1
        load('./results/sampler_probs_norm_c.csv');
    else
    % draw probabilities
    nsamp = 20*10^6;
    sampler = arrayfun(@(prob)gamrnd(prob,1,nsamp,1),capHtmod_c(1,1:(end-1))','UniformOutput',false);
    sampler_probs = cell2mat(sampler');
    sampler_probs_norm_c = sampler_probs./sum(sampler_probs,2);
    % very large file!
    writematrix(sampler_probs_norm_c,'./results/sampler_probs_norm_c.csv');
    end

    CE_GLTHI_rea = zeros(1,7);
    CE_GLTHI_abi = zeros(1,7);
    CE_GHHW_rea = zeros(1,7);
    CE_GHHW_abi = zeros(1,7);
    CE_ST_abi= zeros(1,7);
    CE_ST_rea = zeros(1,7);

    % draw welfare result from each degenerate simplex
    for i=1:7
     file = sprintf('./results/state%d.mat',i);
     load(file);
     CE_GLTHI_rea(1,i) = simulation_output.maintable.CE_GLTHI(1);
     CE_GLTHI_abi(1,i) = simulation_output.maintable.CE_GLTHI(2);
     CE_GHHW_rea(1,i) = simulation_output.maintable.CE_GHHW(1);
     CE_GHHW_abi(1,i) = simulation_output.maintable.CE_GHHW(2);
     CE_ST_rea(1,i) = simulation_output.maintable.CE_ST(1);
     CE_ST_abi(1,i) = simulation_output.maintable.CE_ST(2);
    end
  
    gamma = model_parameters.gamma;
    expu_GLTHI_rea = sampler_probs_norm_c*(arrayfun(@(x)u_cara(gamma,x),CE_GLTHI_rea))';
    CE_GLTHI_rea_mix = arrayfun(@(x)invG_cara(x,gamma),expu_GLTHI_rea);
    
    expu_GHHW_rea = sampler_probs_norm_c*(arrayfun(@(x)u_cara(gamma,x),CE_GHHW_rea))';
    CE_GHHW_rea_mix = arrayfun(@(x)invG_cara(x,gamma),expu_GHHW_rea);
    
    expu_GLTHI_abi = sampler_probs_norm_c*(arrayfun(@(x)u_cara(gamma,x),CE_GLTHI_abi))';
    CE_GLTHI_abi_mix = arrayfun(@(x)invG_cara(x,gamma),expu_GLTHI_abi);
    
    expu_GHHW_abi = sampler_probs_norm_c*(arrayfun(@(x)u_cara(gamma,x),CE_GHHW_abi))';
    CE_GHHW_abi_mix = arrayfun(@(x)invG_cara(x,gamma),expu_GHHW_abi);
    
    expu_ST_abi = sampler_probs_norm_c*(arrayfun(@(x)u_cara(gamma,x),CE_ST_abi))';
    CE_ST_abi_mix = arrayfun(@(x)invG_cara(x,gamma),expu_ST_abi);
    
    expu_ST_rea = sampler_probs_norm_c*(arrayfun(@(x)u_cara(gamma,x),CE_ST_rea))';
    CE_ST_rea_mix = arrayfun(@(x)invG_cara(x,gamma),expu_ST_rea);

    cost = sampler_probs_norm_c*model_parameters.MU(1:end-1,:,1);
    
    results = table(cost, CE_GLTHI_abi_mix, CE_GLTHI_rea_mix, CE_GHHW_abi_mix, CE_GHHW_rea_mix, CE_ST_abi_mix, CE_ST_rea_mix);

    writetable(results,'./results/results_probs_7_c.csv');











